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ABSTRACT 

Particle methods are a ubiquitous tool for solving the Vlasov-Poisson equation in comoving coor¬ 
dinates, which is used to model the gravitational evolution of dark matter in an expanding universe. 
However, these methods are known to produce poor results on idealized test problems, particularly 
at late times, after the particle trajectories have crossed. To investigate this, we have performed a 
series of one- and two-dimensional “Zel’dovich Pancake” calculations using the popular Particle-in- 
Cell (PIC) method. We find that PIC can indeed converge on these problems provided the following 
modifications are made. The first modification is to regularize the singular initial distribution function 
by introducing a small but finite artificial velocity dispersion. This process is analogous to artificial 
viscosity in compressible gas dynamics, and, as with artificial viscosity, the amount of regularization 
can be tailored so that its effect outside of a well-defined region - in this case, the high-density caustics 
- is small. The second modification is the introduction of a particle remapping procedure that peri¬ 
odically re-expresses the dark matter distribution function using a new set of particles. We describe 
a remapping algorithm that is third-order accurate and adaptive in phase space. This procedure pre¬ 
vents the accumulation of numerical errors in integrating the particle trajectories from growing large 
enough to significantly degrade the solution. Once both of these changes are made, PIC converges at 
second order on the Zel’dovich Pancake problem, even at late times, after many caustics have formed. 
Eurthermore, the resulting scheme does not suffer from the unphysical, small-scale “clumping” phe¬ 
nomenon known to occur on the Pancake problem when the perturbation wave vector is not aligned 
with one of the Cartesian coordinate axes. 

Subject headings: cosmology: theory — methods: numerical — dark matter 


I. INTRODUCTION 

Cosmological dark matter (DM) simulations follow the 
evolution of a self-gravitating, collisionless fluid in a co¬ 
ordinate system that expands along with the universe. 
This fluid can either be treated in isolation or coupled to 
a second baryonic fluid that is subject to its own set of 
physical processes. Such simulations cover a wide range 
of length scales, from cosmic-scale (^ 10 Gpc) calcula¬ 
tions that cov er a representative volume of the observable 


universe (e.g. Alimi et al. 20I2|), down to galactic scale 
(^ 10 kpc) calculations that probe the substructure of 
the DM fluid within Milky-Way sized (^ 10^^ ^ 0 ) halos 
(e.g. |Kuhlen et al.||2008|). The results of these simula¬ 
tions aremsed~ToguiHe7interpret, or calibrate virtually 
all current and planned observational efforts to study the 
nature of dark matter and dark energy, including direct 
and indirect dark matter detec tion, galaxy redshift sur¬ 
veys, and weak leasing studies (jKuhlen et al.||2012|). 

The evolution of the DM fluid is governed by the 
Vlasov-Poisson equation in comoving coordinates. While 
it is possible to solve this set of equations using grid - 
based methods in phase space (Yoshikawa et al.||2013), 
this approach is not common in practice due to the 
computational expense of working in high-dimensional 
spaces. Typically, particle methods are used instead. 
These methods are based on a Lagrangian description 
of the Vlasov-Poisson system, in which the initial dark 
matter distribution function is sampled by a finite set of 
particles. The particle positions and velocities are then 


advanced in time using the appropriate equations of mo¬ 
tion derived from the Vlasov-Poisson equation, combined 
with some sort of scheme for computing the forces at the 
particle positions. This force calculation can be carried 
out in a variety of ways. In PIC methods, for example, 
an Eulerian representation of the density is constructed 
from the particle positions by deposition and the Poisson 
equation for the gravitational potential is solved on the 
resulting mesh points. An example of a cosmologic al DM 
code that uses pu re PIC for its force solve is MC^ (Heit- 
mann et al.| | 2QQ5[); other pop ular choices are tree codes , 


e.g. pdkgrav2 ptadel 2001) and 2HOT (Warren 
adaptive PIC codes, e^ ART ( KravtsoV et ah ' 


RAMSES (|Teyssier|2002 ) and Nyx (jAlmgren et al. 2013 ), 
or co mbinations of the two, e.g. GAD GETS ( [Springel 



2005|) and HACC (|Habib et al.||20I2|). While the de¬ 
tails of the force solves difter, all of these codes share the 
underlying particle discretization of the Vlasov-Poisson 
system. Such calc ulations have now b een run with hun- 


dreds of million s (Ali mi et al.||20I2|) to o ver a trillion 
(Ishiyama et al.||2012[ 


Skillman et al. 


20I4[ ) particles. 


T'he accuracy of cosmological observations has ad¬ 
vanced tremendously over the past few decades, and with 
the advent of next-generation instruments like the Large 
Synoptic Survey Telescope in sight, the prospect of 
measuring the basic cosmological parameters of the uni¬ 
verse to_bofcter_thunM%a£CUTa£y has become conceivable 
(e.g. Heitmann et al.||2009 2010|. However, interpreting 
these observations will require theoretical predictions of 
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similar or better accuracy, and it is well-known that stan¬ 
dard particle techniques produce poor results on some 
idealized test problems. A particularly striking exam¬ 
ple is the off - axis plane wav e collapse problem of |Melott 


et ah (1997). Melott et al. set up a plane wave pertur- 


bation that was arbitrarily aligned relative to the Carte¬ 
sian coordinate axes and traced its evolution into the 
non-linear regime using several different particle-based 
numerical schemes. All of the schemes failed to converge 
to the correct solution when sparse particle counts were 
used. Instead of respecting the symmetry of the problem 
setup, non-physical clumping was observed in directions 
perpendicular to the axis of the perturbation. This result 


has since been reproduced by sever al researchers (Heit- 


mann et al. 12005 Hahn et 'aL]|2013 ) 


I'he inability of cosmological dark matter simulations 
to converge on this problem is troubling. It is not clear 
to what extent these problems manifest themselves in 
more realistic calculations. Code comparisons involving 
more realistic initial conditions generally show encour¬ 
aging agreement o n physical observables for dark-matter 


only problems (e.g.|He itmann et al.|2QQ5|[Heitmann et al. 
2QQ8[ Kim et al.||2QM ~ However, all of tnese codes solve 
the Vlasov-Foisson system using particle methods, dif¬ 
fering mainly in the way the force solve is carried out. 
Thus, these comparisons cannot rule out the hypothe¬ 
sis that common assumptions are giving rise to common 
errors. A related issue is the small-scale artificial frag¬ 
ment at ion_obseryedJi]__paxti^e simulations of warm dark 
matter (Wang & White||2007). This fragmentation may 
be related to the “clumping” effect observed in the off- 
axis pancake problem, but this needs to be d emonstrated 
explic itly. Additionally, as pointed out by [Hahn et al. 


(2013), we cannot be sure that cold dark matter simu¬ 


lations do not suffer from the same problems, but are 
simply hidden by the presence of initial power on small 
scales. 

In light of this, a particle method that can be rigor¬ 
ously shown to converge - even on a simple test problem 
like the Zel’dovich Pancake - is desirable. The problem 
of improving the accuracy of particle meth ods on cos¬ 
molog y calculations has been considered by | Hahn et al.| 
(2013), who presented an improved PIC scheme that ex¬ 
ploits the continuity of the infinitesimally thin DM sheet 
in phase space. We take an alternative approach, based 
on th e conv ergence theory for PIC methods from |Wan^ 
et al. (2011 hereafter WMC2011), which focused on the 
application of PIC to electrostatic plasmas. That theory 
shows that PIC methods can converge on the Vlasov- 
Poisson problem provided that 1) the initial particle 
spacing is less than or equal to the cell spacing of the 
Poisson mesh, and 2) a phase-space remapping procedure 
is applied to the particles that prevents the accumulation 
of numerical error in the particle trajectories from over¬ 
whelming the solution. However, this theory assumes 
that the initial distribution function is finite, and hence 
does not apply to the cold dark matter version of the 
Vlasov-Poisson problem, where the initial distribution 
function is usually taken to be a delta-function in veloc¬ 
ity space. Thus, in order to use this theory to guide our 
investigation of the cold dark matter problem, we will 
regularize the initial conditions so that the convergence 
theory in |WMC2QII| applies. A natural way to do so is 
to introduce a finite, artificial velocity dispersion cji to 


the initial DM distribution function. Once solutions to 
these artificially warm versions of the cold dark matter 
problem are obtained, they can used to be used to ob¬ 
tain solutions to the perfectly cold problem by taking the 
limit ^ 0. 

An instructive analogy is to the use of an artificial vis¬ 
cosity in shock-capturing schemes for compressible gas 
dynamics. For Reynolds numbers characteristic of super¬ 
sonic fiows in astrophysics, the molecular mean free path 
of the gas is tiny compared to the other length scales of 
interest, and thus the effects of viscosity should be neg¬ 
ligible. However, neglecting the viscous terms entirely 
leads to numerical instabilities that can contaminate the 
solution far from any shocks. A common solution is to 
introduce an artificially large numerical viscosity that 
smooths out the solution near a shock front over a few 
grid cells. Away from the shock, the integrity of the 
solution is preserved. Likewise, in structure formation, 
the present-day microscopic velocity dispersion of dark 
matter particles is negligibly tiny compared to the bulk 
velocities that arise from gravitational collapse. How¬ 
ever, ignoring the velocity dispersion completely leads to 
resolution-dependent results at the caustic locations, and 
could potentially cause other problems as well. Using an 
artificially large initial velocity dispersion smooths out 
the dark matter caustics over a length scale that is nu¬ 
merically resolvable, preserving the qualitative features 
of the dark matter density distribution without affecting 
the solution outside of the caustics. 

The second mod ification suggested by the error anal¬ 
ysis of WMC2QII is to periodically remap the DM par¬ 
ticles in phase space. This addresses a downside of par¬ 
ticle methods, which is that numerical errors made in 
integrating the particle trajectories tend to compound 
with each successive time step. In electrostatic PIC, 
this phenomenon manifests itself as an exponentially 
growing term in the stability error for the electric field 
(WMC20II Equation 3.12). This term eventually over- 
whelms the solution as the simulation evolves. To miti¬ 
gate this, the distribution function must be periodically 
re-expressed using a new set of particles before integra¬ 
tion errors have a chance to accumulate. This procedure 
is crucial for obtaining quality solutions to plasma prob¬ 
lems for long time integrations. Note that this procedure 
requires the above regularization of the initial data, in 
that we must be able to generate a well-defined distri¬ 
bution function in phase space corresponding to a given 
particle distribution. 

In this paper, we show that this approach improves the 
accuracy of PIC on the Zel’dovich Pancake test problem. 
The structure of the paper is as follows. We begin by 
reviewing the cosmological Vlasov-Poisson system a^ 
the standard PIC algorithm for solving it in Section^ 
In Section we perform a series of convergence stu( 
ies that test PIC’s ability to track the gravitational col¬ 
lapse of a single, plane-wave perturbation to late times. 
We perform two versions of this test, first with only one 
spatial dimension, and second in two spatial dimensions, 
with the axis of the perturbation misaligned with the 
Cartesian coordinate axes. We confirm that standard 
PIC with singular initial data converges poorly at late 
times in ID, and that it suffers from unphysical clump¬ 
ing on the 2D, oblique problem, regardless of the number 
of particles per Poisson cell employed. Next, in Section 
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we consider the same numerical method, using initial 
data that has been regularized via the introduction of a 
finite artificial velocity dispersion. We show that, given 
finite initial conditions, the standard PIC technique does 
converge at 2nd order in ID. We also show that, in the 
limit that the artificial velocity dispersion goes to zero, 
we recover the solution to the problem with perfectly 
cold initial conditions. However, we find that the in¬ 
troduction of regularization does not, by itself, solve the 
problem with artificial clumping on the 2D, oblique prob¬ 
lem. Finally, in Section we introduce particle remap¬ 
ping. We perform this remapping on an adaptive set 
of grids that automatically tracks the evolution of the 
dark matter distribution function in phase space as the 
universe expands. Our results using both regularization 
and remapping are shown in Section!^ We show that, 
when remapping is employed, 2nd order convergence re¬ 
sults are still obtained in ID, and the solution is not 
affected by the artificial clumping present in other cases. 
We summarize our conclusions in Section [3 


2. EQUATIONS AND ALGORITHMS 
2.1. The Vlasov-Poisson system 

We work with a non-dimensional form of the Vlasov- 
Poisson equation in comoving coordinates: 


dt 


^ At + 

a dx 



V + -Vcj) 
a 


dv 


(1) 


Here, f{x,v) is the dark matter distribution function in 
2D-dimensional phase space (x^v) G x where D 
is the number of spatial dimensions under consideration. 
X and V are the us ual comoving position and peculiar 
velocity coordinates (Peebles|1993), related to the proper 
position coordinate r by 


X = a{t) 

V = a{t)x. (2) 

The expansion of the universe is described by the time- 
dependent scale factor a(t). In general, the time evolu¬ 
tion of a{t) is determined by the assumed cosmological 
parameters. In this paper, we specialize to a flat, critical- 
density, matter-only universe, so that 


a{t) = 



(3) 


The gravitational potential (j) obeys the Poisson equa¬ 
tion: 

VV=|;(.-1). (4) 

where p is the comoving dark matter density. We com¬ 
plete the specification of the problem by adopting peri¬ 
odic boundary conditions on the domain x G (0,1)^ and 
supplying appropriate initial conditions, f{x,v^t = 


2.2. Partiele Diseretization 

In particle methods for numerically solving Equation 
[3 the initial distribution function is discretized by a set 
of Lagrangian particles, P. Given /(cc, v, tini), we must 


choose a set of particles with initial positions veloci¬ 
ties Vp, and masses rrip such that 

f{x, V,tini) « T, - K) ■ 

peF 

From these initial coordinates, the particles evolve ac¬ 
cording to the appropriate equations of motion. The tra¬ 
jectory {xp{t), of oach particle can be obtained by 

substituting Equation into Equation [3 The result is 
the following system ofODEs: 

drrip 
dt 

dxp 

dt 

dvp 

dt 

where Qp is the acceleration on particle p induced by 
the dark matter distribution. It is important to note 
that these particles have nothing to do with the actual 
dark matter particles whose statistics are described by 
the Vlasov-Poisson equation; they are simply a set of 
interpolating points from which the distribution function 
can be recovered. 


= 0 


= -V. 


d 1 

- Vp + -Qp, 

n ^ n t' 


(6) 


2.3. The PIC Update 

Once the initial particles have been generated, we ad¬ 
vance their positions and velocities using a standard 
Particle-in-Cell technique (Hockney & Eastwood 1981). 
We write the particle comoving positions and peculiar 
velocities at time t'^ as Xp and where n is the time 
step number. We also assume that we know the gravita¬ 
tional field at the particle positions at the same time 
index, either by performing an initial Poisson solve or 
by storing the result from the previous time step. Given 
that information, we advance the particles to time 
as follows: 

Particle Kick. For our time integration scheme, we 
use the secon d-order and sympleptic K ick-Drift-Kick se¬ 
quence from 
scheme is sel 


Miniati & Colella (2007). Note that this 
starting in the sense that it does not re¬ 
quire any information from time steps prior to n to com¬ 
plete the update - an important feature given our use 
of remapping below. We begin by updating the particle 
velocities to the half-time index n + 1/2: 


V 


tT/T 1 /2 _ 


1 


^n+l/2"P ^n+l/2^P 2 


(7) 


Here, At is the time step, and and are the 

expansion factors computed at times and 
Particle Drift. Next, the particle positions are ad¬ 
vanced using the half-time velocity: 


^n+l_ I 

d^p dyp T 


,n+l/2 P 


t,n+l/2A^. 


(8) 


The final update to the particle velocities cannot be com¬ 
pleted until we compute the new forces at time index 
n -L 1. 

Particle Deposition. To do so, we compute the den¬ 
sity at the mesh points Xi = {i P 1/2) Ax, where i G 
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Fig. 1.— The overall control flow of a PIC calculation 


are the cell indices: 




E l 


Wcic 


Xj — X 


n+1 


Ax 


(9) 


Here, V = Ax^ is the volume of cell i and Wcic{x) 
is the I)-dimensional cloud-in-cell (CIC) interpolating 
function: 

D 

w CIC {x) = Wcic {Xd) , (10) 

d=l 


Wc\g{x) 


1 — |a;|, 0 < |a;| < 1, 
0 otherwise. 


(11) 


Note that in general, we do not use the same mesh spac¬ 
ing for the particle discretization and the Poisson mesh, 
i.e. Ax 7 ^ hx. 

Poisson Solve. The next step is to solve the Poisson 
equation for the gravitational potential at the same grid 
points on which the density is defined. We discretize 


the Laplacian operator using the standard 2D+1 point 
centered difference approximation: 


D 


-E 

d=l 






Ax" 


2a”+i 


- 1) ■ 


( 12 ) 

To solve the resulting linear sy stem, we use the geo metric 
multigrid package in Chombo (Adams et al.|2015), using 
Gauss-Seidel with Red-Black ordering as the smioother. 
We set the solver tolerance to 10“^ - small enough its 
precise value does not affect our convergence results be¬ 
low. After iterating to convergence, the gravitational 
field g = —V(j) can be computed at the same grid points 


= - 


/n+l 


^n+1 


2Ax 


(13) 


Force Interpolation. Next, we interpolate the field 
back to the particle positions using the same interpolat¬ 
ing function as in the deposition step: 


9 ;^" = ^9iViW CIC 


^. _ ^n+1 

Ax 


(14) 


Particle Kick. Finally, we complete the update to 
the particle positions as: 


^n+l ^ 

Up 


^n-\-l/2 


•.n+l P 


,n+l/2 


+ 


1 


a 


n+l 


(15) 


2.4. Time stepping strategy 

In this paper, the size of the time step At is controlled 
by two factors. First, we limit the factor by which the 
background can expand in a single time step. The time 
step associated with the expansion factor is: 


^^exp Gexp 


(i)' 


(16) 


where Cexp is a tunable parameter. 

Second, we enforce a Courant-Friedrichs-Lewy (CFL) 
-type constraint, which limits the distance that a particle 
can travel in a single time step to be some fraction of Ax. 
The time step associated with this constraint is: 


^^part — 


Ax 


part 


maxdvpl)' 


(17) 


In this paper, we take Cpart =0.5. We use the smaller 
of these two time steps to advance the simulation: 


At = min(Atexp, Atpart)- 


(18) 


This ap proach to time s t eppin g is quite similar to that 
in, e.g., Almgren et al. (2013). In the runs presented 
in this paper, the expansion of the background typically 
controls At at early times, while the particle CFL con¬ 
straint controls it at late times. 

Given an initial set of particles, the above equations 
describe how to advance the system to an arbitrary later 
time. The overall c ontr ol flow of the PIC procedure is il¬ 
lustrated in Figure [ 2 ^ We have implemented the above 
algorithm, as well as the remappiM procedure described 
in Section using the Chombcj^ software framework 


https: //commons, lbl.gov/display/chombo 
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for partial differential equations, using standard message 
passing with block-structured domain decomposition for 
parallelization. This code, as well as our data analysis 
scripts, are available onlin^ 


3. THE ZEL’DOVICH PANCAKE 

We now study the convergence of the PIC method 
for plane-wave initial conditions. The gravitational col¬ 
lapse of a plane-wave perturbation in an expanding back¬ 
ground, sometimes called the “Zekdovich Pancake” prob¬ 
lem because of the flattened structures it prod uces, is 
a common test case for cosmological codes (e.g. Bryan 


et al. 119951 [Kravtsov et al.|1997[ Miniati fc Colella|2(JU7[ 

Hopkins |2014|) for several reasons. At early times (before 

any matter parcels cross) an exact solution exists in ID, 
making the problem valuable for code validation. After 
the first shell crossing, the problem still provides a valu¬ 
able test on a code’s ability to deal with strong density 
contrasts and poorly resolved features. Finally, while the 
problem is clearly idealized, it still physically relevant, as 
the collapse of a single Fourier mode forms the basis for 
more complex structure formation calculations. 


3.1. Initial conditions 

The initial distribution function is taken to be an in¬ 
finitesimally thin sheet in phase space: 

f V, tini) — Pz (^5 ^ini) d (^V V ^ (^X ^ ^ini)) • 

The initial density pz {x^t[^i) and velocity Vz {x^ti^i) 
are then computed using the Zel’dovich approximation. 
For example, if the matter parcels of the unperturbed 
DM fluid are labelled by their comoving Lagrangian co¬ 
ordinates g, and we apply a sinusoidal perturbation of 
the form 

S{q) =—sm{k ■ q)k, (20) 

where A is the amplitude and k the wavenumber, then 
(in our adopted Einstein-de Sitter cosmology) the comov¬ 
ing position and peculiar velocity of the fluid at time t 
are 


x{t) = q aAsm{k • q)k 

( 21 ) 

v{t) = adAsm{k ■ q)k. 

( 22 ) 

The corresponding density is 


1 

^ ^ 1 + aA\k\ cos{k • q) ’ 

(23) 

which can be converted to Eulerian coordinates x using 
Equation 

The above equations are valid until the first caustic 
forms, which happens at the expansion factor 

1 

UcausticHj — A|^|* 

(24) 


At that point, the velocity becomes multi-valued and 
the density diverges at the positions corresponding to 
cos{k q) = 1. Prior to that time, we can use the above 
equation to set initial conditions. In particular, we rep¬ 
resent the above initial conditions with a collection of Np 

^ https://bitbucket.org/atmyers/cosmologicalpic 


TABLE 1 

Summary of parameters 

FOR THE SINGULAR, ID RUNS 


Ncells 

Ux 

C^exp 

256 

128 

1.0 X 10-2 

512 

128 

5.0 X 10-3 

1024 

128 

2.5 X 10-3 

2048 

128 

1.25 X 10-2 


equal-mass particles labelled by the Lagrangian coordi- 

n 3 Fog 

qp = {i + 1/2) h,, (25) 

where hx = l/Np is the initial particle spacing (which is 
not necessarily the same as the mesh spacing for the Pois¬ 
son solve. Ax), i G and rup = l/Np is the particle 
mass. We start the integration at time tint, correspond¬ 
ing to expansion factor aini- At that time, the Eulerian 
positions and velocities of the particles are 


^1 = <lp+ «ini-4sin(fc • qp)k 

Vp = ainidini^sin(fc ■ qp)k. 


(26) 

(27) 


Note that after applying a perturbation, neither the 
particle positions nor the velocities are laid out on a 
Cartesian grid. This procedure, in which the Zel’dovich 
approximation is used to construct initial particle posi¬ 
tions and velocities and the initial distribution function 
is perfectly cold, is standard in cosmo logical dark matter 
simulations (e.g. Hahn fc Abel|[ 2 MT ). 


3.2. ID results 

We begin by performing a resolution study on a ID 
problem setup. This is the easiest test for the method 
to pass, since some numerical instabilities may only be 
possible in multi-dimensional problems. However, this 
test will help us verify that we have implemented the 
above scheme correctly, and it may also provide clues 
about the behavior of the higher dimensional case. 

Eor these runs, we set the perturbation wavenumber 
to 27r, the fundamental mode of the computational box. 
We initialize the problem at aini = 1/200 and choose 
the perturbation amplitude so that shell crossing occurs 
at acaustic = 0.1. We then evolve the simulations until 
^stop = I 7 using the scheme described in Section 2.3 The 
simulation results are dumped out every Aa = 0 . 01 , i.e. 
at a = 0.01, 0 . 02 , ..., 0.99, 1.00. We fix the number of 
particles per Poisson cell, where = Np/NceWs^ to 
128, and vary the number of Poisson cells in the prob¬ 
lem domain Aceiis = 1 /Ax over the range {256, 512, 
1024, 2048}. Every time we decrease the mesh spac¬ 
ing, we also decrease the parameter Cexp by a factor of 
2. These parameters are summarized in Table Note 
that our use of 128 particles per cell is something of a 
best-case scenario; most calculations in cosmology use 
significantly fewer particles. The resulting solutions are 
shown at selected times in Eigure As expected, at 
a{t) = 0.1 the particle trajectories cross each other and 
a single dark matter caustic forms. As the particles con¬ 
tinue to evolve in phase space, secondary and tertiary 
caustics form. Each caustic is marked by a peak in the 
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Fig. 2.— The time evolution of the solution to the ID, singular problem. The left panel shows the density, the middle panel the 
gravitational field, and the right panel the potential. The plotted curves are from the run with A^cells = 2048. For the gravitational field 
and the potential, the A/^cells = 1024 run would be indistinguishable from the plotted solutions at the scales shown. The solutions for the 
density, on the other hand, are not converged; see Figure!^ The different lines correspond to the results at different expansion factors: 
solid blue line - a{t) = 0.1; dashed green line - a{t) = 0.3; dashed-dotted red line - a(t) = 0.5; dotted cyan line - a{t) = 0.7. 



X 


Fig. 3.— A zoomed-in view of one of the outer density caustics 
in the ID, singular problem, taken at a(t) = 1. The different colors 
correspond to different resolutions: solid blue line - A^cells = 256; 
dashed green line - A^ceiis = 512; dashed-dotted red line - A^cells = 
1024; dotted cyan line - A^cells = 2048. 


density profile and cusp in the gravitational field. The 
density peaks in the dark matter caustics are clearly not 
converged, as is to be expected with singular initial data 
(see Figure]^. Visually, however, the gravitational field 
and the potential appear to approach well-defined solu¬ 
tions. 

To verify this, we compute the convergence rates for 
the density, gravitational field, and potential. While an 
analytic solution to the Zekdovich Pancake problem ex¬ 
ists prior to acaustic? we are also interested in the con¬ 
vergence behavior at late times, and it is convenient to 
have an error metric that applies to both a < acaustic 
and a > acaustic- Additionally, as we will later want to 
apply our metric to regularized and remapped runs, we 
cannot compare particle quantities directly. Instead, we 
must focus on the grid-defined quantities used in the PIC 
update: p*, and 0*. We therefore use Richardson ex¬ 

trapolation to compute our error estimates. If is one 
of these grid quantities, and is the same quantity 
computed with all the discrete elements (Ax, At, 
Cexp) coarsened by a factor of 2, then we define the rel¬ 


ative solution error as 


= IQ'* - Q 


2h\ 


(28) 


To compare solutions with different numbers of grid 
points, we average the finer solution down to the resolu¬ 
tion of the coarse solution. The order q of the method is 
then 

’ =(tw) ' 

where ||x|| is one of the Li, I/ 2 , or L^o norms of x. In 
what follows, we use the symbols g^, and g^ to de¬ 
note the order of convergence in the density, gravitational 
field, and potential, respectively. 

The resulting convergence rates are displayed as a func¬ 
tion of time in Figure Formally, our adopted PIC 
scheme should be second-order accurate in space, and 
it indeed achieves this up until first shell crossing, at 
least for the gravitational field and the potential. Af¬ 
ter acaustic: howcvcr, thc convergence rates for pi and 
drop dramatically, in all of the norms we consider. The 
solution for pi does not converge at all in the max norm, 
while the solution for converges only slowly. 

In a sense, it is not surprising that we see poor conver¬ 
gence in these metrics after acaustic- With perfectly cold 
initial conditions, there is nothing to limit the peak den¬ 
sities in the dark matter caustics that form after acaustic 
other than the finite numerical resolution employed. In¬ 
deed, we see in Figured that the Loo-norm of p actually 
diverges with g^ < 0 after acaustic- These results are con¬ 
sistent with the analytic theory of caustic structure (in 
the no n-gravitating case) from Shandarin & Zeldovich 
(1989), which states that, after pancake formation, there 


should be a p{x) (x x singularity in the density, and 
therefore a g{x) (x cusp in the gravitational field. 
While this poor convergence in ID is limited to the caus¬ 
tic positions, it points towards more significant problems 
in 2D or 3D calculations, where the symmetry constraints 
of the initial conditions are not automatically met. 


3.3. 2D Results 

In this section, we examine a 2D analog of the non 
axis-aligned p l ane w ave collapse problem considered by 
Melott et al. (1997). This setup, in which the pertur- 
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Fig. 4.— Convergence results for the singular version of the ID pancake problem. The plotted curves show the Richardson-estimated 
order versus the expansion factor a(t). For the solid blue line, the order was estimated using the runs with A^cells = 256, 512, and 1024, 
while the dashed green line shows the same quantity for resolutions A^cells = 512, 1024, and 2048. The top row shows the rates computed 
for the density, the middle row the gravitational field, and the bottom row the potential. Each column shows the rate using a different 
error norm; from left to right, we show the Li, L 2 , and Loo norms, respectively. The expansion factor at which the first caustic forms is 
indicated with a dashed vertical line. 


TABLE 2 

Summary of parameters 

FOR THE SINGULAR, 2D RUNS 


-^cells 

nx 

Dexp 

128^ 

1/2^ 

2.0 X 10-2 

128^ 

12 

2.0 X 10-2 

128^ 

162 

2.0 X 10-2 

2562 

1/2^ 

1.0 X 10-2 

2562 

12 

1.0 X 10-2 

2562 

162 

1.0 X 10-2 


bation wave vector is not aligned with any of the coor¬ 
dinate axes, is known to produce unphysical clumping 
along directions in which the density should be constant 
given the symmetry of the initial conditions. In this sec¬ 
tion, we confirm that this effect is also present in 2D 
using our standard PIC method. We set the perturba¬ 
tion wavenumber k = {kx^ky) to be (2, 5) in units of the 
fundamental mode. Otherwise, the initial conditions are 


the same as in our ID resolution study: Uini = 1/200, 
^icaustic = 0.1, and Ustop = 1- Throughout the rest of this 
paper, we refer to these initial conditions as the “oblique” 
pancake. 

We perform six runs in total, as summarized in Table 
We vary the number of cells in the Poisson mesh 
over {128^, 256^}, and run the problem with Ux = 1/2^, 
1, and 16^ particles per Poisson cell. The particles are 
initially arranged within the Poisson cells in the usual 
way. In the run with 1/2^ particles per cell, we coarsen 
the Poisson mesh by a factor of two and place particles 
at centers of each cell in the resulting coarsened mesh. 
In the run with 16^ particles per cell, we refine the mesh, 
instead. 

The resulting particle x and y positions at a = 1 are 
shown in Figure We only display the high resolution 
results in Figure^ the low resolution results are qual¬ 
itatively similar. By the symmetry of the problem, the 
particle density should be constant along directions per¬ 
pendicular to the axis of the perturbation. However, all 
of the runs show signs of small-scale fragmentation that 
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do not obey this symmetry requirement, and thus cannot 
be valid solutions to Vlasov-Poisson for the given initial 
conditions. This effect is not limited to runs with sparse 
particle counts, and is even present in the run with the 
16^ particles per Poisson cell - a much larger number 
than typically used in realistic dark matter simulations. 
We conclude that simply using very large particle counts 
with the basic PIC scheme is not sufficient to prevent 
these errors. 

This clumping effect is present in the grid-deposited 
density field, as well. In Figure we construct a 2D 
density field by depositing the particles onto a 256^ mesh 
using cloud-in-cell deposition (Equation]^. This Figure 
was constructed using the run with 16^ particles per cell; 
the other runs produce an even stronger effect. Figure 
shows the corresponding density field from the 16^ par¬ 
ticle per cell, high resolution run, where we generate the 
density by deposition onto a 512^ mesh. Both figures 
shows signs of unphysical fragmentation along the dark 
matter caustics. 

It is possible that this artificial fragmentation is related 
to the divergence of the d ensit y at the caustic locations 
in ID observed in Section [3^ To test this hypothesis, 
we introduce a regularization procedure that limits the 
maximum density at the caustic positions. The results 
of this test are described in the next section. 


4. THE REGULARIZED PANCAKE 
4.1. Initial Conditions 

The above initial conditions are singular in the sense 
that the initial distribution function is a delta function in 
velocity space, centered at v = Vz- With delta function 
initial conditions, the solution to the Zel’dovich pancake 
problem at late times contains singularities in the den¬ 
sity, and corresponding discontinuities in the derivative 
of the gravitational field, at the positions of the dark mat¬ 
ter caustics. When solving this problem numerically with 
PIC or other particle methods, however, the peak densi¬ 
ties obtained in the numerical solution will be determined 
by the finite amount of resolution employed. In effect, 
the mass deposition kernel, which spreads out the mass 
in a given particle over a few grid cells, provides a reg¬ 
ularization that removes these singularities; however, it 
does so in a relatively uncontrolled way that is explicitly 
resolution-dependent. As a consequence, with perfectly 
cold initial data, it is not possible to separate numeri¬ 
cal errors, which should improve with mesh refinement, 
from parameters of the numerical model, which should 
have well-defined effects on the solution once a converged 
answer is found. 

Instead of letting the mesh spacing provide the reg¬ 
ularization, we investigate an alternative regularization 
procedure that removes the singularities in the solution 
by introducing a finite initial velocity dispersion, to 
equation We choose a Gaussian form for the velocity 
profile, so that the initial distribution becomes: 


f{x,V,tini) = 


27rai 


D/2 


exp 


(v V ^ini)) 


2a,' 


ip) 


1 + aA|/g| cos{k • q) 


(30) 


This process makes the initial conditions more closely 
resemble those used in electrostatic PIC calculations, for 
which the expected convergence rates can be demon¬ 
strated. In the limit a, ^ 0, Equation is recovered. 
Using these modified initial conditions, we then evolve 
the regularized version of the Zel’dovich Pancake prob¬ 
lem using the same PIC method as before. To generate 
the initial particles, we sample Equation with a set 
of particles that are initially laid out on a cell-centered, 
Cartesian grid in phase space. We label this grid Uq. 
The computational domain stretches from 0 to 1 in each 
physical dimension and — U to U in each velocity dimen¬ 
sion. The domain bounds in velocity space are chosen 
so that f{x,v) is negligibly small outside of the problem 
domain. The number of points in Qq in the position and 
velocity dimensions are Nx and Ny, respectively. The 
locations of the cell centers are thus cc, = {i l/2)hx 
and Vj = {j + l/2)hy — U, where {i,j) G (Z^,Z^), 
and {hx^hy) = {l/Nx,2V/Ny) are the cell spacings in 
position and velocity space. We place one particle in 
each cell (i, j). For a given initial distribution function, 
/(cc, v,tini), the initial discretization is completed by as¬ 
signing each particle p G P a mass 


'mp = f{xl,vl,tini)h^h^. (31) 

where Xp and Vp are the initial position and velocity of 
particle p. For computational efficiency, we discard par¬ 
ticles with masses less than 10“^^; experimentation re¬ 
veals that tightening this value does not meaningfully af¬ 
fect our convergence results. We illustrate the difference 
between our “regularized” and “singular” initial parti¬ 
cle layouts in Figure We also show, in Figure a 
sample spectrum of the particle masses obtained through 
this procedure. Specifically, these were the actual parti¬ 
cle masses used in the a, = 1, A'ceiis = 256 calculation 
described in Table Note that there is long tail of low- 
mass particles that contribute little to the overall distri¬ 
bution function - this is a consequence of our selecting 
a relatively low particle mass threshold and a relatively 
large domain in velocity space. It is likely that, by re¬ 
laxing these parameters, significantly fewer particles can 
be used without degrading the solution. 

As an aside, although our regularized initial conditions 
are “warm” in a sense, they should not be confused with 
“Warm Dark Matter”. Our initial velocity dispersion is 
artificially large and introduced for numerical purposes 
only. It has nothing to do with the physical velocity 
dispersion of a hypothetical DM candidate. Simulations 
of structure formation in Warm DM universes in fact use 
“singular” initial conditions by our definition, but they 
suppress power in fluctuations below some free-streaming 
scale associated with the rest mass of the assumed dark 
matter particle. 


4.2. ID Results 

To begin, we fix the regularization parameter at a, = 1. 
We otherwise repeat the calculation in Section us¬ 
ing the same values for aini, acaustic, ^stop, k, and Cexp- 
For our lowest resolution run, we use A^ceiis = 256, 
Nx = 512 = Ny = 512. T hese choices were motivated by 
the convergence theory of WMC2QTTj which states that, 
for convergence. Ax <= hx and hx = 0{hy) (in our di¬ 
mensionless formulation of the Vlasov-Poisson problem). 
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Fig. 5. — Particle x and y positions at a{t) = 1.0 on the singular version of the oblique pancake problem. All runs use a 256^ Poisson 
mesh. The top row shows a run with a sparse particle count of 1 particle every 4 Poisson cells. The middle row uses 1 particle per cell, 
and the bottom uses a very high particle count of 256 particles per cell. For all three runs, the leftmost panels show the entire problem 
domain, and the view zooms in as you move to the right. All runs show signs of small-scale fragmentation that is inconsistent with the 
symmetry of the problem. 


While this doesn’t require that we find it 

convenient to enforce this condition. This amounts to 
employing a relatively fine base grid, leading us to over¬ 
resolve the problem in velocity space, but this is not much 
of an issue in this 1+lD problem. We also perform runs 
with the resolution increased to A^ceiis = 512, 1024, 2048, 
and 4096, and the other elements of the discretization 
refined accordingly. That is, when we refine the Pois¬ 
son mesh, we also refine Ny, and decrease the time 
step parameter Cexp, as before. These run parameters 
are summarized in Table The resulting solutions are 
shown in Figure The Sect of the (in this case quite 


large) artificial velocity dispersion to smooth out the den¬ 
sity peaks over some comoving length scale, so that the 
peak densities in the caustics now have well-defined max¬ 
imum values. As shown in Figure 10, the density caustics 
now appear to be converging withincreased mesh reso¬ 
lution. 

We also examine the convergence rates for the den¬ 
sity, gravitational field, and potential on the regularized 
version of t he problem. The rates are computed as in 
Section |3.2| and shown in Figure Overall, the con¬ 
vergence behavior is much better than on the singular 
version of the problem. The density order is quite noisy 
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Fig. 6.— The density field at a(t) = 1 from the low resolution, 
computed by deposition onto a 256^ mesh. 


singular version of the oblique pancake problem. The density has been 


TABLE 3 

Summary of parameters for the regularized, ID 

PANGAKE RUNS 


CTi 

a 

-^cells 

iVa, 

Ny ^ 

Dexp 

{1,1/2,. 

.,1/32} 

64 

128 

128 

4.0 X 10-2 

(1,1/2,. 

.,1/32} 

128 

256 

256 

2.0 X 10-2 

(1,1/2,. 

.,1/32} 

256 

512 

512 

1.0 X 10-2 

(1,1/2,. 

., 1/32} 

512 

1024 

1024 

5.0 X 10-3 

(1,1/2,. 

., 1/32} 

1024 

2048 

2048 

2.5 X 10-3 

{1,1/2,. 

.,1/32} 

2048 

4096 

4096 

1.25 X 10-2 


^ The notation {1,1/2,..., 1/32} means we have varied 
the parameter cr^ over the stated range. 

^ The number of cells in the initial particle mesh in the 
X direction. 

The number of cells in the initial particle mesh in the 
V direction. 

at early times, but the rates for the potential and gravi¬ 


tational field are close to 2 for all the times we consider. 
The convergence rate for the density is slower than 2 at 
late times, but the density is not divergent as it was on 
the singular problem. 

4.3. The Double Limit 

In Section |4.2[ we showed that for a given value of 
cTi, our PIC method converges at late times on the ID, 
regularized pancake problem. Of course, in dark matter 
simulations, we are really interested in the case where 
(Ji = 0. Fortunately, we can exploit the fact that we 
are able to obtain converged solutions to the regular¬ 
ized problems to also obtain a solution to the singular 
problem by considering the limit as ^ 0. This proce¬ 
dure gives us a way of separating numerical errors, which 
should improve with mesh refinement if our method con¬ 
verges, from regularization error, i.e. the fact that our 
caustics are artificially smoothed out by the artificial ini¬ 
tial velocity dispersion. In this section, we illustrate this 
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Fig. 7.— The density field at a{t) = 1 from the high resolution, singular version of the oblique pancake problem, using the same color 
scale as Figure The density has been computed by deposition onto a 512^ mesh. 
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Fig. 8.— A schematic demonstration of the difference between regularized and singular initial conditions. Left panel - the initial phase- 
space discretization for our regularized runs, for Aa:/hx = 2. The red crosses show the initial particle locations, and the black dots mark 
the centers of the Poisson mesh. The particle masses are computed by sampling the initial distribution function. Right - the initial particle 
positions for our singular runs, also for = 2. The particles all have equal masses, and have been displaced from their initial positions 

using the ZePdovich approximation. 
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10 “^^ 10 “^° 10 “^ 10 “"^ 


Fig. 9.— A histogram of the particle ma sses obtained through 
the sampling procedure described in Section [4.1| We have used 50 
logarithmically-spaced mass bins, and show /p, the fraction of the 
total number of particles in each bin. To generate this plot, we 
have used CTj = 1.0, V = 6.0, and take n Nx = Ny = 512. All other 
parameters are the same as in Section |3.2| 



Fig. 10.— A zoomed-in view of one of the outer density caustics 
from the ID, regularized problem with cr^ = 1, taken at a(t) = 1. 
The different colors correspond to different resolutions: solid blue 
line - A^cells = 256; dashed green line - A^cells = 512; dashed-dotted 
red line - A^cells = 1024; dotted cyan line - N^eWs = 2048. 


process on our ID problem setup. We vary the quantity 
ai over the range {1, 1/2, 1/4, 1/8, 1/16, 1/32}. For 
each value of cr^, we increase the mesh resolution until 
we obtain a converged solution (see Table for a full 
list of the run parameters). Next, we consioCT the limit 
of the converged solutions to the regularized problems 
as (Ji ^ 0. This concept of a “double limit” - letting 
bo th the mesh sp acing and go to zero - was inspired 
by Krasny (1986), who used a similar procedure to re¬ 
move singularities in vortex methods for incompressible 
fluid flows. The results are shown in Figure These 
figures show that the regularized solutions approach the 
singular solution as —> 0. The main difference between 
the singular and regularized runs is that the peak caus¬ 
tic densities are limited in the regularized runs, and the 
structures in the inner caustics are smoothed out over 


some comoving length scale. 

In particular, as is decreased, we resolve more and 
more of the internal structure of the dark matter density 
distribution. To make this more quantitative, we asso¬ 
ciate a comoving length scale d with each solution as 
follows. For each value of <j^, we compute the comoving 
distance (at a{t) = 1) between the positions of the outer¬ 
most caustics that are present in the perfectly cold run, 
but not present in the corresponding regularized run. To 
identity the peak positions, we use a standard clump find¬ 
ing algorithm. The results are as follows. For the largest 
value of = 1, even the outermost caustic positions 
differ slightly from the results obtained in the perfectly 
cold limit, while the internal caustics have been com¬ 
pletely washed out, such that none of the caustics with 
positions within d ~ 0.032 of the center in the cold prob¬ 
lem are present in the warm problem. For = 1/4, the 
outermost caustic positions match precisely, while caus¬ 
tics within d ~ 0.008 of the center in the cold problem 
have smoothed out. For = 1/16, the corresponding 
length scale is d 0.0013. 

Outside of the caustics, the effect of on the solutions 
is small. This finding suggests an alternative approach 
to perfectly cold initial conditions in DM simulations. 
Instead of using singular initial conditions, which make 
code validation difficult, one instead decides on a length 
scale below which one will not believe the answer, and 
then chooses the corresponding initial artificial velocity 
dispersion. Given that velocity dispersion, one can rig¬ 
orously demonstrate convergence at the desired order of 
accuracy. 

To make the concept of the double limit more quanti¬ 
tative, we define the following error metric, analogous to 
the Richardson extrapolated error introduced in Section 

Mi 

ey<T) = ||5?(<T)-gbA2)||. (32) 

That is, for each value of the mesh spacing h, we compute 
the norm of the difference between the solution with the 
regularization parameter = a and the one with ai = 
a/2. We use the L 2 norm to compute the error metric, 
and focus on the solution for the gravitational field. The 
results for e^{a) for all the runs in Tableware shown in 
Figure[^ for the same expansion factors snown in Figure 

There are two conclusions to be drawn from this plot. 
First, it shows that the quantity e^{a) approaches a well- 
defined value for each pair of velocity dispersions. By 
^ceiis = 2048, e^(cr) has leveled off for all of the values of 
ai considered. Second, it confirms that the regularized 
solutions themselves converge as ^ 0. 

There are several conclusions that can be drawn from 
this exercise. First, it shows that our method of regular¬ 
izing the initial conditions through an artificial Gaussian 
velocity dispersion does not affect the dark matter struc¬ 
ture outside the caustic positions, provided ai is chosen 
small enough. Second, it gives us reason to believe that 
the solution to the singular problem obtained via the 
standard PIG method is in fact the correct answer to the 
ID pancake problem after Ucaustic- Because the density 
field diverges in the perfectly cold limit, the standard 
method of validating the solution through examining the 
convergence rate is not available. By considering the 
limit of a series of regularized problems, we can have con¬ 
fidence that the solution obtained to the singular prob- 















The Convergence of Cosmological PIC Schemes 


13 




Fig. 11.— The time evolution of the solution to the ai = 1 version of the ID, regularized problem. The left panel shows the density, 
the middle panel the gravitational field, and the right panel the potential. The plotted curves are from the run with A^cells = 2048; the 
A^cells = 1024 run would be indistinguishable from the plotted solutions at the scales shown. The different lines correspond to the results 
at different expansion factors: solid blue line - a(t) = 0.1; dashed green line - a(t) = 0.3; dashed-dotted red line - a(t) = 0.5; dotted cyan 
line - a(t) = 0.7. 
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Fig. 12.— The same as Figure^ but for the = 1 version of the regularized problem. The blue line compares the result for the 
-^cells = 256, 512, and 1024, while the green line shows the same quantity for resolutions A^cells = 612,1024, and 2048. 
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TABLE 4 

Summary of parameters for the 

REGULARIZED, 2D PANCAKE RUNS 


CTi 

-^cells 

Nx ^ 

iv„ b 

Dexp 

1/16 

128^ 

256 

64 

2.0 X 10-2 

1/16 

256^ 

512 

128 

1.0 X 10-2 


lem is correct. Third, since we are not able to solve the 
oblique version of this problem directly using cold initial 
conditions, regularization may provide a way to obtain 
approximate solutions to the singular problem in more 
than one spatial dimension. We test this hypothesis in 
the next section. 


4.4. 2D Results 

We now compute the solution to a reg ular ized version 
of the 2-dimensional problem in Section [T3l We do two 
runs in this section, which are summarized in Table 
Both runs have set = 1/16, and construct Hq suoi 
that Nx = 256, and Ny = 128. In this case, the total 
number of cells in the phase space domain on which the 
particles are initialized is Vj x N^. The number of cells 
in the Poisson mesh is fixed at Vceiis = 128^. All of the 
other parameters are the same as in Section 3.3[ The 
resulting particle positions for the higher resolution are 
shown in Fi gure [l5l We also display deposited fields, 
as in Section [33| for both the low (Figure [T^ and high 
(Figure [T^ resolutions runs. 

Unfortunately, despite the improved convergence rates 
in ID, the regularized problem suffers from the same 
clumping problem as the singular problem for the oblique 
pancake setup. However, there is an additional refine¬ 
ment we can introduce to the basic PIC procedure. In 
the next section, we introduce particle remapping, and 
show that, in concert with the regularization of the ini¬ 
tial conditions described in this section, it can improve 
the performance of PIC on the oblique problem. 


5. REMAPPING 

In the error analysis of electrostatic PIC in |WMC2QII[ 
the stability error for the electric field contains a term 
that grows exponentially with time. Given long enough 
integration times, this error term can become large 
enough to significantly degrade the quality of the nu¬ 
merical solution. A common strategy for controlling this 
errors is to periodically restart the problem with a new 
set of particles, before the accumulated errors in the par¬ 
ticle trajectories become too large. Such “regridding” or 
“remapping” techniques have been applied successfully 
to particle schemes for fluid dynamics, such as vortex 


(iVadlamani et al. 20041 Koumoutsakos 1997 Cottet & 

Koumoutsakos||2000 

Chaniotis et al.||2002p, 

and to FKJ 

in the context of pJ 

asma physics 

(IDenavit 

1972 

Chen 

et al.||2008 Wang et 

al.||20il 2012 

)■ 


distribution function back onto a regular mesh in phase 
space, initialize a new set of particles, and restart the 
calculation. In principle, the remapping mesh can uni¬ 
form. However, repeatedly re-sampling the distribution 
function on a uniform mesh leads to a large number of rel¬ 


atively weak particles in the tails of the distribution func¬ 
tion - a poor use of the available resolution elements. Fur¬ 
thermore, in the cosmological context, remapping serves 
an additional role. We have regularized the initial con¬ 
ditions by introducing a finite initial velocity dispersion, 
ai. However, as the scale factor a{t) increases, a{t) will 
scale as (Ji/a{t). Thus, if we employ the above proce¬ 
dure an a uniform mesh, the artificial velocity spread in 
the distribution function will be resolved more and more 
poorly as time advances, undoing the natural adaptivity 
of Lagrangian schemes. 

To remedy these issues, we remap on an adaptive hier¬ 
archy of meshes that automatically adjusts the velocity- 
space resolution on the finest level so that the distri¬ 
bution is resolved by approximately the same number of 
particles throughout the simulation. The AMR hierarchy 
is defined as a set of cell-centered, Cartesian meshes in 
phase space. We label the meshes where 0 < ^ < ^max 
are the level numbers. The coarsest mesh, Uq, covers the 
entire phase-space domain and is the same as the mesh 
used for generating the initial particles. The finer meshes 
are unions cell-centered rectangles and in general cover 
only a portion of the problem domain. The mesh spac¬ 
ing of level ^ is related to the mesh spacing on the next 
coarsest level by hi = hi-xjr^Qi^ where Vref is the AMR 
refinement ratio. The valid region of an AMR level is 
defined as the region not covered by finer levels, and the 
composite grid is defined as the union of the valid regions 
of each level. 

We describe the remapping algorithm in detail below, 
first using a single, uniform mesh to represent /(cc,i;), 
second using a fixed hierarchy of levels, and finally us¬ 
ing the full, adaptive hierarchy. Note that, as presently 
constructed, the remapping procedure requires the reg¬ 
ularization of the initial conditions - there must be a 
well-defined distribution function to generate the parti¬ 
cles from. 


5.1. Remapping on a Uniform Mesh 

As an illustration, consider remapping a set of parti¬ 
cles on Qq only. With only one lev el, our remapping al¬ 
gorithm is identical to |WMC2QII[ The remapping step 
proceeds as follows: 

Deposition. From the particle data, we represent 
f{x,v) on Uq by deposition. The position and ve¬ 
locity coordinates of cell centers on Uq are the same 
as those used for initialization, Xi = {i ^ l/2)hx and 
Vj = {j + lf2)hy — U, with («,j) G (Z^,Z^). The val¬ 
ues of / on the mesh points are thus: 


/«= E (^) 


Xi — x^ 

hx 


Wa 


'^3 - 

hy 


(33) 

where F = is the phase-space cell volume and 

Wii x) is the 2D-dimensional versio n of the thi r d-ord er 
accurate interpolating function from Monaghan (1985): 


D 


Wi (x) = 11 Wi (Xd) 


(34) 


d=l 


1 - 


Wi{x) = 


2,^, . ikp, 0<N<1, 

i (2 - (1 - |a:;|), 1 < ia:;i < 2, (35) 

0 otherwise. 
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Fig. 13.— The effect of letting cr^ —)■ 0 on the obtained solutions. The solid lines show the converged density (top row), gravitational 
field (middle row), and potential (bottom row) from our regularized runs at a(t) = 1. The black dotted lines are taken from a very high 
(-^cells = 16, 384) resolution singular calculation. The regularized results are taken from the runs with A^cells = 4096 - the highest resolution 
available. The solutions are sufficiently converged that the A^cells = 2048 solutions would be indistinguishable from the plotted curves on 
this plot. The left column (blue) shows CTj = 1, the middle column (green) shows CTj = 1/4, and the right column (red) shows a = 1/16. 


High-order interpolation is necessary beca use one orde r 
of accuracy is lost during the remap step (WMC2011). 
Thus, for the overall scheme to have second-order accu- 
racy, a third-order or higher interpolating function must 
be used. 

Positivity Preservation. A consequence of high- 
order interpolation is that / is not guaranteed to be pos¬ 
itive for all cells in Qq. To account for this, a correction, 
A/, must be applied to the distribution function after it 
is deposited onto the mesh. We accomplish this using a 
mass transfer procedure that redistributes matter to cells 


with / < 0 from neighboring cells in proportion to each 
neighbor’s current value of /. We write the undershoot 
in cell {ij) as dfij: 

5fij (36) 

and initialize A/ = 0. For every cell with negative Sf, we 
compute the capacity of each neighboring cell Ci+rn,j+n 
as 

Ci+m,j+n = max(0,, (3^) 

where the notation {i + m, j + n) refers to any neigh¬ 
boring cell in the redistribution zone. We perform this 
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Fig. 14.— Our error metric e^{a) computed for all runs in Table 1^ The blue curve shows e^(cr = 1) - that is, the L 2 norm of the change 
in the solution for the gravitational field when a changes from 1.0 to 0.5, as a function of h. The green curve corresponds to cr = 0.5, the 
red curve a = 0.25, cyan is cr = 0.125, and magenta is 0.0625. The values in this plot were computed at a(t) = 0.1, (top left), a{t) = 0.3, 
(top right), a(t) = 0.1, (bottom left), and a(t) = 0.1, (bottom right). 



Fig. 15.— Particle x and y positions at a{t) = 1 for the ai = 1/16 version of the regularized oblique pancake problem. The scales shown 
are the same as in Figure 
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Fig. 16.— The solution for the density at a, = 1 from the low-resolution calculation of the oblique pancake problem using regularized 
initial conditions with cr^ = 1/16. The color scale is the same as in Figure]^ Particle remapping was not turned on for this run. For details 
of the parameter choices we adopted, see text 
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Fig. 17.— The solution for the density at a = 1 from the high-resolution calculation of the oblique pancake problem using regularized 
initial conditions with Gi = 1/16. The color scale is the same as in Figure]^ Particle remapping was not turned on for this run. For details 
of the parameter choices we adopted, see text. 


redistribution over a window of M cells in each direc¬ 
tion. The value of M must be chosen so that it matches 
the number of cells involved in the particle interpolation 
stencil: in this case, A/" = 4. The total capacity is 

neighbors 

^tot ~ ^ ^ ^i+m’,j+n’' 

The correction to / within the redistribution window 
is incremented as 

^fi+mj+n += — ^fij‘ (39) 

Ftot 

This procedure is repeated for each non-positive cell in 
and the correction is then applied to /. Note that 
there is no guarantee that / will be positive after just one 
pass of this algorithm. In general, we must iterate until 
/ is positive everywhere. The full positivity preservation 


algorithm looks like: 
while max(/) > 0 do 
Set A/ = 0 
for all (i,j) G do 

if fij < 0 then 

for all neighbors {i + m^j + n) do 

^fi+mj+n += — ^fij- (40) 

Ftot 

end for 
end if 
end for 

jnew ^ yold ( 4 ;^) 

end while 
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Fig. 18.— The overall control flow of our PIC algorithm, with 
remapping. 

In practice, we find that 2 or 3 iterations is usually suf¬ 
ficient to ensure that the positivity of the distribution 
function is preserved. 

Particle generation. Finally, we generate a new set 
of particles by sampling the remapped distribution at 
every cell (i, j) in Oq and initializing a particle with mass 

mp = fijh^h^. (42) 

As in the problem initialization, we discard particles with 
masses less than 10“^^. 


deposition. However, we find that this technique leads to 
spurious oscillatory features in the distribution function 
when used in our remapping algorithm (see Figure p^. 
Instead, we find it necessary for each particle to remem¬ 
ber the spacing of the mesh on which it was generated. 
We label these spacings and e^. 

With that caveat, the remapping step proceeds as fol¬ 
lows: 

Deposition. We deposit the particles in each set onto 
the appropriate level. For example, on 



for (i, j) G f2o, (43) 

where the expression is evaluated for every cell (i^j) in 
f2o, not just the valid cells. The quantities 6x^ and Sv^ 
are: 

6x^ = max(ha,rfgf, e^) 

Sv^ = max(h^rfgf, e^). (44) 

That is, we use whichever is larger, the mesh spacing 
of the target mesh or the mesh spacing with which the 
particles were generated. We find that this interpolation 
scheme results in a smooth representation of the under¬ 
ling distribution function even when adaptive refinement 
is employed. 

The corresponding expression for level 1 is: 



for {i,j) e (45) 


The second term in this expression accounts for the 
particles in Pq that deposited some of their mass on the 
region covered by fti. The indices are eoarsened 

versions of (i, j): 


5.2. Remapping on a Fixed Hierarehy 

Next, consider the case with a hierarchy of levels, {f2o, 
..., We impose the constraint that these 

levels be properly nested by a buffer of ribuff = 4 fine-level 
cells; this choice is determined by the above interpolation 
stencil. For now, we assume that the levels are fixed in 
time and known in advance; we consider the full AMR 
case in the next section. We must now compute / on the 
composite grid. To do so, we first partition the particles 
P into sets {Pq, Pi, ..., P^max-i} according to what level 
they will be deposited on. We perform this partition by 
associating each particle with the finest level that con¬ 
tains its entire interpolation stencil; particles that lie too 
close to a coarse/fine boundary for their entire clouds to 
be contained on the fine level go on the next coarser level. 

Remapping on a hierarchy of meshes raises the addi¬ 
tional issue of what mesh spacing to use in the deposi¬ 
tion step: the spacing of the mesh that the particles are 
being interpolated to, or the spacing of the mesh the par¬ 
ticles were generated from. It is common for AMR PIC 
schemes to use the spacing of the target mesh for density 


^ = 


3 = 


^0 

F ref. 

jo 

F ref 


'^D-l 
F ref 

jP-l 

F ref 


(46) 


Note that by adopting a 4-cell proper nesting require¬ 
ment, we avoid the case where this coarse/fine correction 
can involve more than 2 levels. This process continues 
up to ^max- In general: 




5v^ 


, u-i 

I JiCjC^ 


for (ij) G 


(47) 


The result of applying this procedure over levels 
..., ^ representation of / that fully conserves 

matter over the composite grid. 

Positivity Preservation. The positivity preserva¬ 
tion step is similar to the above. The difference is that 
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Fig. 19.— A ID demonstration of our method for depositing particle quantities on refined meshes. In both panels, an initial sinusoidal 
function f(x) (the black curve) has been sampled by a set of 32 particles, such that hx = 1/32. These particles were then deposited onto 
a 64 cell mesh (Ax = 1/64). This mesh was used to generate a new set of particles. In the left panel, the coarser spacing hx was used for 
deposition, while in the right panel, the fine cell spacing Ax was used instead. The coarse spacing results in a much smoother representation 
of the underlying function. 


when applying positivity preservation algorithm to level 
we must sometimes draw mass from cells that lie on 
level £ -\- 1 and £ — 1. For convenience, we accomplish 
by augmenting each level with a later of of N ghost cells. 
When we apply the correction to /^, we also apply it to 
the subsets of levels £-\-l and £—1 that are covered by the 
level £ ghosts. To coarsen the correction, we use simple 
arithmetic averaging, and to refine it, we use piecewise- 
constant interpolation. 

Particle Generation. The particle generation step 
proceeds as in the single-level case, except that we gener¬ 
ate one particle for every valid cell in the hierarchy. That 
is, we do: 

for £ = 0, ..., £top - 1 do 

for all {i,j) e valid do 

Initialize particle with mass nip = ffjh^h^■ 

end for 
end for 


5.3. Remapping with AMR 

In practice, the AMR hierarchy is not known at the 
beginning of a remapping step. Instead, the AMR level 
structure must be self-consistently built up from the par¬ 
ticle distribution and a suitable set of refinement criteria. 
In this paper, we use the condition that a cell is tagged 
for refinement whenever the value of ffj exceeds the pre¬ 
determined threshold /thresh = 0.1. Finer levels are then 
generated in such a way that they cover these tags, and 
are consistent with our proper nesting requirement. The 
other input to the process is the maximum number of al¬ 
lowed levels, ^max- We then build up the AMR represen¬ 
tation of / using the standard bootstrapping procedure: 

for £top = 0, ..., ^max - 1 do 
Partition particles into {Pq, 

Deposit particles on {f2o,..., ^4op} 
if "^top ^max 1 then 

Regrid: 

for £ = 0, ..., £top do 

Tag cells on where > /thresh 

end for 


for £ = 1, ..., £top + 1 do 
Generate new 

end for 
end if 
end for 

The positivity preservation and particle generation 
steps are the same as for the fixed hierarchy case. The 
result of this procedure is an AMR representation of /. 
See Figure for an illustration in the 1 + ID case. 
Note that, since the particle remap is a purely local pro¬ 
cess, it is not necessary (nor would it be advisable in 3 
-f 3D) to store the full phase space distribution at once. 
However, in the present work we confine our attention 
to low-dimensional problems, so we store the full / for 
simplicity and for illustrative purposes. 


5.4. Refinement Criteria 


The goal of our use of AMR is to resolve the artifi¬ 
cial velocity dispersion by a roughly constant number 
of particles during the remap step, even as decreases 
with the inverse of the expansion factor. We achieve this 
by fixing the velocity-space domain boundaries and the 
cell spacing of the base grid and adding additional lev¬ 
els of refinement to the remapping mesh as the universe 
expands. We compute the number of AMR levels by re¬ 
quiring that cr(a) is always resolved by at least N^- cells. 
If the cell spacing of the base grid in velocity space is 
then we perform the remapping on an AMR hierarchy 
with 


log {N^hy/cJi) 
log (rref) 


(48) 


levels of refinement, where Tref is the AMR refinement 
ratio. 


5.5. Summary 

We show the updated control flow of PIC with remap¬ 
ping in Figure After initialization, we b egin PIC 
time integration using the time step in section 2.4 We 
continue until we reach either the final expansion fac¬ 
tor, Ustop, oi* the expansion factor associated with the 
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Fig. 20.— An example of a 1 + ID distribution function represented on an adaptive hierarchy of grids. The level 0 grids have been 
omitted for clarity, while the level 1, 2, and 3 grids have been over-plotted in white, gray, and black, respectively. 


TABLE 5 

Summary of parameters for the regularized and 

REMAPPED ID PANGAKE RUNS 


CTi 

-^cells 

Na 

Adremap 

Ax 

Ny 

Uexp 

1.0 

256 

2 

0.01 

128 

128 

1.0 X 10-2 

1.0 

512 

4 

0.01 

256 

256 

5.0 X 10-3 

1.0 

1024 

8 

0.01 

512 

512 

2.5 X 10-3 

1.0 

2048 

16 

0.01 

1024 

1024 

1.25 X 10-3 


first remap. The remapping can be applied either in 
fixed interval in a or every fixed number of time steps. 
For each remap step, we compute the finest allowed level 
^max using Equation 48, and perform the remap on 
..., }• This process continues until Ustop is reached. 


6. THE REGULARIZED, REMAPPED PANCAKE 

6.1. ID Results 


TABLE 6 

Einal energy error as a fungtion of resolution 


Acells 

Singular ^ 

Regularized ^ 

Remapped ^ 

256 

6.3 X 10-4 

3.2 X 10-4 

1.2 X 10-5 

512 

1.5 X 10-4 

8.2 X 10-5 

1.5 X 10-4 

1024 

4.2 X 10-5 

2.0 X 10-5 

7.1 X 10-6 

2048 

1.2 X 10-5 

6.1 X 10-6 

3.3 X 10-6 


^ The displayed quantity is the energy conservation 
error e, evaluated at a(t) = 1. For the regularized and 
remapped columns, we have used the runs with ai = 1. 


In this section, we examine the performance of the 
PIC scheme with remapping on the ID pancak e pro blem 
with regularized initial conditions (see Section [iTT] ) . Ba¬ 
sic PIC converged on this problem (provided the initial 
conditions were regularized) so we must verify that we 
can also obtain second order convergence with remap¬ 
ping enabled. For this run, we choose = 1, and we 
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arrange our remapping meshes so that this velocity dis¬ 
persion is always resolved by = 8 cells. The initial 
phase-space grid is N^xNy^ where = l/h^ = 512, 
Ny = 1/hy = 512, and the Poisson mesh is Vceiis = 256 
cells across. By a{t) = 1, the mesh spacing in velocity 
space has been refined by a factor of 2^ on the finest level. 
We apply the remapping every Aaremap = 0-01 (i.e. at 
a = 0.01, 0.02, etc... ) and dump out f {x,v), p{x), 9 {x)^ 
and (j){x) after every remap. We show the resulting so¬ 
lution at representative times in Figures through 
In addition to the density, gravitational fidd, and poten¬ 
tial, we have also displayed the full 1 + ID phase-space 
distribution function at the same expansion factors. As 
expected, most of the phase space cells are empty as thus 
do not need to be represented by a particle. Visually, the 
regularized and remapped solutions appear quite similar 
to the regularization-only solutions, with the effect of 
being to smooth out the density at the caustic locations 
so that rigorous convergence becomes possible. 

To investigate the convergence rate, we run an anal¬ 
ogous set of problems where we vary all of hy^ 

and At over the same values as in Section 14.21 Addi¬ 
tionally, we also increase No- by a factor of two as the 
resolution increases. The resulting convergence rates are 
plotted in Figure The convergence rates at late times 
are much improved compared to the singular runs. Ad¬ 
ditionally, the convergence rates for the density are im¬ 
proved beyond the regularization -only case. The noisy 
convergence rates seen in Section |4.2| at early times for 
the density are gone. The reason for this behavior is that 
the remapping procedure constrains the density field to 
be sampled by several particles per Poisson cell every¬ 
where, which results in a smoother representation of the 
density. Additionally, the convergence rates for the den¬ 
sity are better at late times, as well, with rates close 
to 2 being obtained in all the norms. The faster-than- 
second-order convergence rates in the potential are a con¬ 
sequence of the positivity preservation procedure - they 
are not present if this part of the algorithm is disabled. 

It is important to note that, while the above remapping 
scheme conserves the total / exactly (modulo our choice 
to discard particles with masses less than 10“^^), it is not 
guaranteed to conserve total energy, due to discretization 
error in performing the remapping on a mesh with finite 
hx and hy. However, if enough resolution is used, these 
errors can be made negligibly small. To demonstrate 
this, we evaluate the overall energy conservation error in 
our runs as a function of time both with and without 
remapping. We define the error in energy conservation 
as follows. If our code conserved energy exactly, then 
the particles w ould obey the Layzer-Irvine equation (e.g., 
Peeblesfl993 ): 

|[a(T + [/)] = -dT, (49) 



X 


Fig. 21.— A zoomed-in view of one of the outer density caustics 
from the ID, regularized and remapped problem with cr^ = 1, taken 
at a(t) = 1. The different colors correspond to different resolutions: 
solid blue line - A^cells = 256; dashed green line - A^cells = 512; 
dashed-dotted red line - A^cells = 1024; dotted cyan line - A^cells = 
2048. 

This quantity would be zero in the case of perfect 
energy conservation. In Figure we compare e as 
a function of time for the singular pancake run, the 
(Ti = 1 regularized run, and the = I regularized and 
remapped run. We use the highest resolution data avail¬ 
able (Aceiis = 2048). Because of our use of a variable 
time step, neither the singular nor the regularized runs 
conserves energy exactly, despite our use of a sympletic 
integrator. However, the overall error in the energy con¬ 
servation is small - at then end of all three runs, we have 
conserved energy to better than 0.002 %. Interestingly, 
the energy error is actually largest at the end of the sin¬ 
gular run, which suffers a spike in the energy error at 
first caustic formation. 

The energy error is practically the same between the 
regularized and remapped runs - an indication that, at 
this resolution, the energy error is dominated by some¬ 
thing other than discretization error during the remap. 
For our lower resolution runs, this is not the case. In 
Table we compare the energy error a{t) = I between 
the singular run, the regularized run with = I, and 
the regularized and remapped run with di = 1 for four 
different resolutions. We see that on the Aceiis = 256 
run, the energy error is actually the largest O.I %) 
in the run with remapping; at this resolution, the dis¬ 
cretization error in the remap step is clearly significant. 
By Aceiis = 512, e in the remapped run has become com¬ 
parable to that in the singular run, and at still higher 
resolution, it becomes comparable to the error in the 
regularized run. 

6.2. 2D Results 


where T = ^ rriivf is the total kinetic energy of the 
particles and U = ^ is their total gravitational 

potential energy. We use this equation to define the en¬ 
ergy conservation error as e: 

_ a{T + U)-ao{To + Uo) + j;^Jda 


Finally, we repeat the oblique pancake problem with 
both regularization and remapping turned on. The run 
parameters are summarized in Tablef^ We set di = 1/16 
and vary Aceiis over {128, 256} . We use velocity-space 
mesh refinement during the remap phase so that di is 
always resolved by at least No- = 1 cells, up until a maxi¬ 
mum refinement level of 2 is reached. We remap in linear 
increments in the expansion factor, with Auremap = O.OI. 
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Fig. 22. — The energy conservation error e (see text for definition) 
as a function of a(t) for the singular run (blue solid line), the regu¬ 
larized run (green dashed line), and the regularized and remapped 
run (red dashed-dotted line) for the ID pancake problem. In this 
figure, we have used the highest resolution data available, for which 
-^cells = 2048. The variation in the ener^ error with respect to 
mesh refinement is summarized in Table Im The vertical dashed 
line marks the expansion factor at which the first caustic forms. 


The resulting particle positions are displayed in Figure 
26 Unlike Figures 1^ and which show the particles 


positions at a = 1 for the singular and regularized but 
not remapped runs, Figure ^E\ shows no sign of particle 
clumping along the caustic axes. Instead, a consequence 
of the remapping procedure is that the particle spacings 
are forced to be relatively constant throughout the popu¬ 
lated regions of phase space. The information about the 
density is contained in the particle masses, instead. Fig¬ 


ures 


28 and shows the grid-deposited density at a = 1 


for 128^ and 256^ run, respectively, where we have con¬ 
structed the 2D density field by depositing the particles 
onto a 256^ and 512^ mesh using cloud-in-cell deposition. 
The unphysical clumping visible in Figures and is 
no longer visible in these plots. 

It is important to note that this improvement is not 
simply a matter of using more particles per Poisson cell. 
Our regularized, rem apped calculation of the oblique 
pancake (Section 6.2) uses approximately 380 particles 
per Poisson cell on average (the precise number varies 
with time), almost identical to the number used in the 
regularized but not remapped run 385), and similar 
to the 256 used in the corresponding singular calcula¬ 
tion. Simply sampling singular initial conditions with 
more particles in physical space does not lead to a con¬ 
vergent method - our results suggest that the initial con¬ 
ditions must be regularized and the extra particles must 
be arranged in the correct way in phase space. 

However, while using many particles per cell is not suf¬ 
ficient to prevent artificial fragmentation, our results do 
suggest that it may be necessary in our scheme to carry 
around more particles per Poisson cell than is currently 
standard practice. In 3D, even using a relatively modest 
number of particles in each velocity space direction (s^, 
4, comparable to our low resolution runs in Section]^, 
would require 4^ = 64 particles where one would exist 
in the comparable singular problem. Clearly, optimizing 
both the basic PIC scheme and our remapping code for 


TABLE 7 

Summary of parameters for the regularized and 

REMAPPED 2D PANGAKE RUNS 


CTi 

-^cells 

Na 

Aciremap 



Cexp 

1/16 

1282 

I 

0.01 

256 

128 

2.0 X 10-2 

1/16 

2562 

1 

0.01 

512 

128 

1.0 X 10-2 


such high particle-per-cell counts will be necessary before 
our approach can be used for large-scale problems. 

We examine our oblique solutions more quantitatively 
in Figure Here, we plot the max norm of g±_ over 
the max norm of where g_\_ is the component of the 
gravitational acceleration perpendicular to k and g\\ is 
the corresponding component parallel to k. To make a 
clean comparison between runs, we use the gravitational 
force as computed on the mesh points by the PIC al¬ 
gorithm in computing this quantity (Equation p^. In 
the exact solution, this quantity should be zero, Wt due 
to unavoidable discretization errors in the representation 
of the initial conditions, this will not be precisely true. 
These errors are largest at the caustic positions, due to 
the large gradients in the density there. Thus, in all three 
runs, the degree of asymmetry in the force increases dra¬ 
matically after Ucaustic- 

To an extent, these errors are an unavoidable conse¬ 
quence of representing a plane-symmetric quantity on a 
Cartesian grid when that plane is misaligned with the co¬ 
ordinate axes. A stable numerical scheme, however, will 
limit the ability of these errors to affect the subsequent 
particle trajectories. As the error analysis of |WMC2Qll| 
shows, however, PIC is not stable to these errors un¬ 
less remapping is applied. To investigate this, we have 
also shown in Figure the growth in the perpendicular 
component of the deposited velocity, u_l, where the de¬ 
posited velocity was computed by CIC deposition as in 
Equation!^ to make a clean comparison between our dif¬ 
ferent runs. As with g±^ we find that v± is the largest at 
the caustic positions. We also find that, in both our low 
resolution and high resolution models, the application of 
the remapping procedure reduces these off-axis velocities 
by almost an order of magnitude. Note that while v± is 
small in absolute terms, both the perpendicular and the 
parallel components of the deposited velocity should go 
to zero at the caustic positions, and thus any error in 
the velocity is relatively large there. Since many parti¬ 
cles spend a large amount of time the caustics (because 
they are turning points of the dark matter distribution 
function), these small off-axis velocities can lead to sig¬ 
nificant clumping along the axis of the perturbation, if 
nothing is done to mitigate them. 

7. CONCLUSIONS AND CAVEATS 

In this paper, we have investigated the convergence of 
PIC schemes on a variety of “Zekdovich Pancake” se¬ 
tups. Our conclusions are as follows. In one spatial di¬ 
mension, the perfectly cold problem does not converge at 
the formal order of the method at late times, due to sin¬ 
gularities in the solution for the density at the positions 
of the dark matter caustics. This is true regardless of 
the number of particles per Poisson cell employed. Once 
these singularities are removed through the introduction 
of a finite initial artificial velocity dispersion, PIC can 
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Fig. 23.— The time evolution of the solution to the ai = 1 version of the ID, regularized and remapped problem. The left panel shows the 
density, the middle panel the gravitational field, and the right panel the potential. The plotted curves are from the run with A^cells = 2048; 
the A^cells = 1024 run would be indistinguishable from the plotted solutions at the scales shown. The different lines correspond to the 
results at different expansion factors: solid blue line - a(t) = 0.1; dashed green line - a(t) = 0.3; dashed-dotted red line - a(t) = 0.5; dotted 
cyan line - a(t) = 0.7. 


converge at the expected order on the ID problem. Ad¬ 
ditionally, by taking the limit as ^ 0, we have shown 
that we approach the solution obtained in the perfectly 
cold limit, except that the peak density are limited by 
the velocity dispersion. This gives us confidence that, in 
ID, the results obtained with singular initial data are in 
fact valid solutions to the Vlasov-Poisson equation in the 
zero-temperature limit. 

On the 2D, oblique pancake problem, however, the 
basic PIC method with cold initial data does not con¬ 
verge to a valid solution of the Vlasov-Poisson equation. 
Instead, the solutions exhibit small-scale fragmentation 
that does not improve with mesh refinement. This prob¬ 
lem is not solved by increasing the number of particles 
per Poisson cell, nor is it solved by using the artifi¬ 
cially warm initial data. However, the particle remap¬ 
ping algorithm we describe does improve the accuracy of 
PIC on this problem, reducing the off-axis component of 
the gravitational acceleration and velocity significantly 
when applied 100 times over between a(t) = 1/200 and 
a(t) = 1. 

Our results show that remapping is primarily impor¬ 
tant at the caustic positions. This is because that is 
where the gravitational field varies most rapidly with 
position - a consequence of the sharply-peaked density 
structure. Thus, errors in the particle trajectories com¬ 
pound particularly quickly at those locations. Addition¬ 
ally, in more than one spatial dimension, you don’t have 
any control over the directionality of those errors, which 
leads to the clumping observed in the non-remapped 
runs. This fact suggests a possible improvement to our 
remapping scheme: instead of applying it globally to the 
entire distribution function, we could apply it only at 
those spatial locations where the density contrast is large. 

Given that, in ID, the limiting case of the regularized 
solutions agrees with the solution obtained using cold ini¬ 
tial data, and given that the introduction of does not, 
by itself, prevent artificial fragmentation in 2D, a possi¬ 
ble conclusion is that the artificial velocity dispersion is 
not necessary in and of itself, but is primarily important 
in that allows the use of a plasma-style remapping algo¬ 
rithm. If this is the case, it may be possible to construct 
a remapping algorithm that does not require regulariza¬ 
tion, by e.g. exploiting the continuity of the distribution 


function in phase space as in Hahn et al. (2013). How¬ 
ever, whether designing such an algorithrn is possible is 
an open question. To use the remapping algorithms that 
currently exist, regularization is necessary. 

There are a number of caveats to our work that bear 
mentioning. The first is that we have restricted our at¬ 
tention to various configurations of the Zel’dovich Pan¬ 
cake problem. We have not yet applied our technique 
to other, more realistic test problems or to a full struc¬ 
ture formation calculation. Second, we have restricted 
ourselves to working in one- and two-dimensional spaces 
only. Work towards creating a three-dimensional version 
of our remapping algorithm in progress. Finally, while 
our scheme does use adaptivity in velocity space, we have 
not yet introduced AMR in the spatial directions, ei¬ 
ther for solving the Poisson equation or for generating 
the particle positions during remapping. Such a modi¬ 
fication would clearly be advantageous, in that it would 
concentrate both the particles and the resolution of the 
force grid towards the regions where high resolution is 
most needed. However, this was not necessary for the 
relatively simple problems considered in this paper. 
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Fig. 24. — Time evolution of the distribution function for the ai = 1 version of the regularized and remapped problem. Top Left - f{x, v) 
at the expansion factor a(t) = 0.1. Top Right - same, at a{t) = 0.3. Bottom Left - same, at a(t) = 0.5. Bottom Right - same, at a(t) = 0.7. 
The bounds of the images in velocity space have been adjusted slightly in each plot. 

, for data analysis and plotting. 


REFERENCES 


(Hunter[2QQ7 


Adams, M., Colella, P., Graves, D. T., et al. 2015, Chombo 
Software Package for AMR Applications - Design Document, 
Tech, rep., Lawrence Berkeley National Laboratory 
LBNL-6616E 

Alimi, J.-M., Bouillot, V., Rasera, Y., et al. 2012, ArXiv e-prints, 
arXiv:1206.2838 

Alimi, J.-M., Bouillot, V., Rasera, Y., et al. 2012, ArXiv e-prints, 
1206, 2838 

Almgren, A. S., Bell, J. B., Lijewski, M. J., Lukic, Z., & Van 
Andel, E. 2013, ApJ, 765, 39 

Bryan, G. L., Norman, M. L., Stone, J. M., Gen, R., & Ostriker, 
J. P. 1995, Computer Physics Communications, 89, 149 

Chaniotis, A. K., Poulikakos, D., & Koumoutsakos, P. 2002, 
Journal of Computational Physics, 182, 67 


Chen, Y., Parker, S. E., Rewoldt, G., et al. 2008, Physics of 
Plasmas, 15, 055905 

Cottet, G.-H., & Koumoutsakos, P. D. 2000, Vortex Methods: 
Theory and Practice (Cambridge University Press) 

Denavit, J. 1972, Journal of Computational Physics, 9, 75 

Habib, S., Morozov, V., Finkel, H., et al. 2012, ArXiv e-prints, 
1211, 4864 

Hahn, O., Sz Abel, T. 2011, Monthly Notices of the Royal 
Astronomical Society, 415, 2101 

Hahn, O., Abel, T., & Kaehler, R. 2013, Monthly Notices of the 
Royal Astronomical Society, 434, 1171 

Heitmann, K., Higdon, D., White, M., et al. 2009, The 
Astrophysical Journal, 705, 156 

Heitmann, K., Ricker, P. M., Warren, M. S., & Habib, S. 2005, 
The Astrophysical Journal Supplement Series, 160, 28 


f{x,v) 











26 


Myers, Colella, & Van Straalen 


Li L2 



Fig. 25.— The same as Figure^ but for the ai = 1 version of the regularized, remapped problem. The blue line compares the result for 
the Na = 2, 4, and 8 runs, while the green line is for No- = 4, 8, and 16. 
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Fig. 27.— Top - A measure of the asymmetry in the gravitational force as a function of time for the singular (dashed blue lines), 
regularized without remapping (dash-dotted green lines) and regularized with remapping (solid red lines) versions of the oblique pancake 
problem. Low resolution results are shown using thick lines, and high resolution results are shown using thin lines. The quantity on the 
y-axis is the maximum off-axis (that is, perpendicular to the perturbation axis) force as computed on the mesh in units of the maximum 
on-axis (along the perturbation) force. Bottom - a similar plot that shows the growth in the off-axis component of the deposited velocity 
field; see text for details. 
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Fig. 28.— The solution for the density at a = 1 from the low-resolution calculation of the oblique pancake problem using both regularized 
initial conditions and remapping. The color scale is the same as in FigureFor details of the parameter choices we adopted, see text. 
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Fig. 29.— The solution for the density at a = 1 from the high-resolution calculation of the oblique pancake problem using both regularized 
initial conditions and remapping. The color scale is the same as in FigureFor details of the parameter choices we adopted, see text. 
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